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Cosmic Ray particle production 

J. Ranft a * 

a INFN, Lab. Naz. del Gran Sasso, 1-67010 Assergi, Italy 

The status of some popular models to simulate hadronic and nuclear interactions at Cosmic Ray energies is 
reviewed. The models predict the rise of all the hadronic and nuclear cross sections with energy and a smooth 
(logarithmic) rise of average multiplicities, rapidity plateaus and average transverse momenta with the energy. 
Big differences are found between model predictions partly already at energies, where collider data are available. 
It is argued, that at the highest energies data of the Cosmic Ray cascade can only be reliably interpreted by 
sampling the cascade using more than one model. The importance is stressed to put more effort into the models 
and especially a better understanding of the minijet component at the highest energies. Likewise, experimental 
data on particle production are needed at the highest possible energies, to guide the models. 
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1. Introduction 

The extension of models for multiparticle pro- 
duction in hadron-hadron, hadron-nucleus and 
nucleus-nucleus collisions to be used for the sim- 
ulation of the Cosmic Ray cascade up to Ei a b = 
10 21 cV (corresponding to «Js = 2000 TeV) is 
needed to prepare for the Auger experiment jl] 
as well as for a reliable interpretation of present 
experiments like Agasa g] and Flys Eye [|| , which 
present data in the EeV energy region. The 
need for careful comparisons of hadron produc- 
tion models was stressed at the International Cos- 
mic Ray Conference in Roma 1995. Following 
this, such a code comparison in the energy re- 
gion of interest to the Kaskade experiment j| 
was presented by members of the Kaskade exper- 
iment From this code comparison it became 
clear, that already in the knee region of the Cos- 
mic Ray energy spectrum important differences 
exist between the models and that these differ- 
ences might change the interpretation of certain 
Cosmic Ray results. Here we will discuss the sta- 
tus of some of these models, discuss the minijet 
component, present typical comparisons to Col- 
lider data, present some characteristics of hadron 
production up to Ei a b — 10 21 eV and finally com- 
pare some results obtained simulating the cosmic 
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ray cascade using different models. 

2. The present status of some event 
generators used for Cosmic Ray cascade 
simulations 

The presently dominant hadron production 
models used for the simulation of the Cosmic Ray 
cascade are constructed on the basis of multi- 
string fragmentation, they use Gribov-Regge and 
Gribov-Glauber theory, to construct the multi- 
string production in hadron-hadron and nuclear 
collisions. Most of the models use minijets as an 
important mechanism for particle production at 
high energies. 

The DPMJET-II event generator based on the 
two-component Dual Parton Model (DPM) was 
described in detail The extension of this 

model up to energies of y/s — 2000 TeV was 
reported this year, the resulting model will be 
refered to as DPMJET-II. 3. The extension is 
done by calculating the minijet component of 
the model using new parton distribution func- 
tions, the GRV-LO parton distributions || and 
the CTEQ4 parton distributions fl(i|| , which are 
both available in a larger Bjorken-x range than 
the MRS(D-) parton distributions, which were 
the default in DPMJET-II. 2. These new par- 
ton distributions describe the structure function 
data measured in the last years at the HERA Col- 
lider. DPMJET-II. 3 descibes well minimum bias 
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hadron and hadron jet production up to present 
collider energies. It is also demonstrated, that 
the model performs as well as the previous one 
DPMJET-II.2 for hadron production in hadron- 
nucleus and nucleus-nucleus colisions. DPMJET 
is used for the simulation of the Cosmic Ray cas- 
cade within the HEMAS-DPM code U used 
mainly for the MACRO experiment Jl2| . 

The SIBYLL model |ll| is a minijet model and 
has been reported to be applicable up to Ei a b = 
10 20 eV. However, the EHQL Q parton struc- 
ture functions used for the calculation of the mini- 
jet component might , after the HERA experi- 
ments, no longer be adequate. It is known, that a 
significant updating of SIBYLL is planned for the 
next year. SIBYLL is the most popular model for 
simulating the Cosmic Ray cascade in the USA. 

VENUS, a very popular model applied origi- 
nally for describing heavy ion experiments, is now 
the leading event generator within the Corsika 
Cosmic Ray cascade code jlj|. VENUS is appi- 
cable up to Ei a \, = 5xl0 16 eV. It has been re- 
ported |l6| , that the introduction of minijets into 
VENUS has been planned, this will allow to apply 
VENUS up to higher energies. 

QGSJET [0 is the most popular Russian 
event generator used for Cosmic Ray simulations. 
It is based on the Quark Gluon String (QGS) 
model, this model is largely equivalent to the 
DPM. QGSJET also contains a minijet compo- 
nent and is reported to be applicable up to Ei a t, 
= 10 20 eV. 

HPDM |l8| is based on parametrizations in- 
spired by the DUAL Parton Model, it is reported 
to be applicable up to Ei a b = 10 20 eV, however- 
some of the parametrizations might become un- 
reliable above Ei a t, = 10 17 eV. HPDM was origi- 
nally used as event generator within the Corsika 
cascade code. 

MOCCA @ is an empirical model employing 
a succesive splitting algorihm. It was reported 
to be applicable up to E^ = 10 20 eV. Since the 
model does not contain minijets, its predictions 
at the upper energy end might differ significantly 
from the other models. 

In Table 1 we present some characteristics of 
the models. The Gribov-Regge theory is applied 
by three of the models. The pomcron intercept 



for SIBYLL is equal to one, SIBYLL is a mini- 
jet model using a critical pomeron, with one soft 
chain pair, all the rise of the cross section results 
from the minijets. In the models with pomeron 
intercept bigger than one, we have also multi- 
ple soft chain pairs, already the soft pomeron 
leads to some rise of the cross sections with en- 
ergy. Minijets are used in three of the models, it 
is believed, that minijets are necessary to reach 
the highest energies. All models contain diffrac- 
tive events. Secondary interactions between all 
produced hadrons and spectators exist only in 
VENUS, DPMJET has only a formation zone 
intranuclear cascade (FZIC) between the pro- 
duced hadrons and the spectators. Only three of 
the models sample properly nucleus-nucleus col- 
lisions, the other two models replace this by the 
superposition model, where the nucleus-nucleus 
collision is replaced by some hadron-nucleus col- 
lisions. The residual projectile (and target) nuclei 
are only given by two of the models. 

Table 1. Characteristics of some popular mo- 
dels for hadron production in Cosmic Ray cas- 
cades. (VEN = VENUS, QGS = QGSJET, SIB 
= SIBYLL, HP = HPDM, DPM = DPMJET) 
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3. The calculation of the minijet compo- 
nent 

The input cross section (before the unitariza- 
tion procedure applied by the models) for semi- 
hard multiparticle production (or minijet pro- 
duction) ah is calculated applying the QCD im- 
proved parton model, the details (for DPMJET) 
are given in Ref. pp|-[25|] . 
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fi{x,Q 2 ) are the structure functions of partons 
with the flavor i and scale Q 2 and the sum i,j 
runs over all possible flavors. To remain in the 
region where perturbation theory is valid, a low 
p± cut-off p± thr is used for the minijet compo- 
nent. Since the HERA measurements, the struc- 
ture functions are known to behave at small x 
like l/x a with a between 1.35 and 1.5. The 
minijet production is dominated by very small x 
values, therefore the minijet cross section calcu- 
lated with the new structure functions rise very 
steeply with energy. We found already 1993 p5| 
with the MRS[D-] structure function |^6| at the 
LHC energy a minijetcross section about 10 times 
larger than the total cross section at this energy. 

Such large minijet cross sections are inconsis- 
tent and wrong: The input minijet cross sec- 
tions ah, which one puts into the unitarization 
scheme are inclusive cross sections normalized to 
Tlminijets&inel , where n m imjets is the multiplicity 
of minijets. The physical processes, which con- 
tribute to this inclusive cross section are 2 — > n 
parton processes. 2 — ► n processes give a con- 
tribution to ah equal to na-i^ n . If one treats 
this huge cross section as ah in the usual way 
in the eikonal unitarization scheme one replaces 
it by n/2 simultaneuos 2^2 parton processes, 
this is the inconsistency. What one should re- 
ally use in the unitarization, but what we do not 
know how to compute reliably at present would 
be ah = J2 n a 2^n- The way to remove this incon- 
sistency is to make in the two component DPM 
the threshold for minijet production p±thr energy 
dependent in such a way, that at no energy and 
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Figure 1. Total, inelastic and elastic pp and pp 
cross sections from DPMJET-II.3 as function of 
the center of mass energy y/s. The model results 
obtained using the GRV-LO parton distributions 
U and the CTEQ4 parton distributions (To) are 
compared to the Donnachie-Landshoff fit for the 
total cross section \M and to data |28|-5o| Ell PI 



for no PDF the resulting ah is much bigger than 
the total cross section. Then at least we have 
a cross section, which is indeed mainly the cross 
section of a 2 — > 2 parton process at this level, the 
parton-parton scattering with the largest trans- 
verse momentum. We can get back to the real 
2 — > n processes and recover the minijets with 
smaller transverse momenta via parton shower- 
ing. One possible form for this energy dependent 
cut off is [E5f : 



p ± thr = 2.5 + 0.12[lg 10 (v / I/\/io)] 3 
[GcV/c], V«o = 50GeV. 
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Figure 2. The inelastic cross section (J p -Air calcu- 
lated by DPMJET-II.3 (as well as the ones from 
VENUS and SIBYLL according to Ref. f§) as 
function of the laboratory collision energy (from 
0.02 TeV up to 1.E9 TeV) compared to experi- 
mental data collected by Mielke et al. p7|. 



The resulting <Jh are smaller or not much larger 
than the total cross sections resulting after the 
unitarization for all PDF's. 

There are further features of the minijet com- 
ponent worth mentioning. One uses as first de- 
scribed in |Q at p±thr the continuity requirement 
for the soft and hard chain end p± distributions. 
Physically, this means, that we use the soft cross 
section to cut the singularity in the minijet p± 
distribution. But note, that this cut moves with 
rising collision energy to higher and higher p± 
values. This procedure has besides cutting the 
singularity more attractive features: 

(i) The model results (at least as long as we do 
not violate the consistency requirement described 
above) become somewhat independent from the 
otherwise arbitrary p± cut-off. This was already 
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Figure 3. Pseudorapidity distributions of charged 
hadrons produced in nondiffractive pp collisions 
at Vs = 0.2, 0.54, 0.9 and 1.8 TeV. The 
DPMJET-II.3 results are compared with data 
from the UA-5 Collaboration |}8| and from the 
CDF Collaboration SI. 



demonstrated with DTUJET90 [|| and cut-offs 
of 2 and 3 GeV/c. 

(ii) The continuity between soft and semihard 
physics is emphasized, there is no basic differ- 
ence between soft and semihard chains besides the 
technical problem, that perturbative QCD allows 
only to calculate the semihard component. 

(iii) With this continuity in mind we feel free 
to call all chain ends, whatever their origin in the 
model, minijets, as soon as their p± exceeds a 
certain value, say 2 GeV/c. 



4. Comparing the models to data at accel- 
erator and collider energies 

Each model has to determine its free param- 
eters. This can be done by a global fit to all 
available data of total, elastic, inelastic, and sin- 
gle diffractive cross sections in the energy range 
from ISR to collider experiments as well as to the 
data on the elastic slopes in this energy range. 
Since there are some differences in the hard par- 
ton distribution functions at small x values re- 
sulting in different hard input cross sections we 
have to perform separate fits for each set of par- 
ton distribution functions. After this stage each 
model predicts the cross sections also outside the 
energy range, where data are available. In Fig. [I] 
we plot for DPMJET II. 3 the fitted cross sections 
obtained with two PDF's together with the data. 
Furthermore we compare the total cross sections 
obtained with the popular Donnachie-Landshoff 
fit [^7|. For applications in Cosmic Ray cascade 
simulations we need in particular the hadron-Air 
cross section, in Fig 
to Mielke et al 



we compare data according 
7 1 with the cross sections ac- 
cording to three models. At low energies all mod- 
els are describing these data rather well. At high 
energies we observe however small differences bet- 
ween the models. 

At higher energies (and in non-single diffractive 
pp collisions) there are pseudorapidity distribu- 
tions from the UA-5 Collaboration ]38[ ] and from 
the CDF Collaboration [|9| . In Fig.||a very good 
agreement is found of DMJET-II.3 with these 
data. Still very often there is and was always 
(see Fig.||) a disagreement of the models with the 
UA-5 data at the highest pseudorapidity values. 
The models predict systematically more particles 
at the largest pseudorapidities of the experiment. 
This disagreement (if the data would be correct) 
would of course be of importance, if one is inter- 
ested in Cosmic Ray cascades, where the particle 
production in the fragmentation region is of main 
interest. Fortunately, a new independent mea- 
surement of the pseudorapidity distribution in the 
collider energy range became available recently 
40 1 . In Fig. |J the comparison with this new data 
is presented and we find a remarkable agreement 
with DPMJET-II.3 in the large pseudorapidity 
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Figure 4. Pseudorapidity distributions of 

charged hadrons produced in nondiffractive pp 
collisions at y/s = 0.63 TeV. The DPMJET-II.3 
results are compared with recent data from Harr 
etal. |§. @. 



region. In FigJ^ we present the comparison (from 
Ref . H ) of multiplicity distributions according to 
5 models with the data from the UA-5 Collabora- 
tion |38j . Most of the models describe at least the 
high multiplicity tail of the data reasonably well, 
however the multiplicity distribution occording to 
the SIBYLL model is everywhere rather far from 
the data. We turn to collisions with nuclei. In 
Fig. g the comparison of DPMJET-II.3 is with 
the rapidity distribution of charged hadrons in 
p-Ar collisions at 200 GeV. In Fig. ^| we com- 
pare with the rapidity distribution of negatively 
charged hadrons in central S-S and S-Ag colli- 
sions. 

At least in models with a minijet component 
we expect good agreement with data on trans- 
verse momentum distributions. In FigJ| we com- 
pare hadron jet production in DPMJET-II.3 with 
data from the CDF-Collaboration Q. The jets 
from the model arc found out of the Monte Carlo 
events using a jet finding algorithm with the same 
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Figure 5. Multiplicity distribution of charged 
hadrons from nondiffractive p — p collisions at y/s 
= 540 GeV. The data are from the UA-5 Col- 
laboration |}8| . The comparison with 5 models is 
from IB. 



Figure 6. Charged particle rapidity distribution 
for p-Ar interactions. The DPMJET-II.3 results 
are compared with data |M2j . 



parameters like the one used by the experiment. 
With a minimum bias Monte Carlo event gener- 
ator it is of course not possible to obtain good 
statistics on the total transverse energy range of 
the experiment. We find good agreement of the 
jets in the model with the data up to E± = 30 
GeV/c. The transverse momentum distribution 
in a large pj_ region was determined by the UA- 
1-MIMI Collaboration 45 1 . In Fig.^J we compare 
DPMJET-II.3 results with the parametrization 
of the data given by this experiment and we find 
a good agreement. 

In Fig.[l0| we compare average transverse mo- 
menta as obtained from DPMJET-II.3, QGSJET 
and SIBYLL as function of the cms energy yfs 
with data collected by the UA-1 Collaboration. 
This plot gives at the same time the DPMJET 
predictions for the average transverse momenta 
up to yfs = 2000 TeV and the predictions of the 
two other models up to y/s = 100 TeV. While at 
energies where data exist all models agree rather 
well with each other and with the data, we find 



completely different extrapolations to higher en- 
ergies. We should note, this are just the three 
models with a minijet component. But it seems, 
that in spite of the minijets the average transverse 
momentum in QGSJET becomes constant at high 
energies, while it continues to rise in DPMJET. 
For me the rise of the average transverse momen- 
tum in DPMJET is connected with the fact, that 
with the new parton structure functions since the 
HERA measurements really the minijets domi- 
nate very much all of hadron production at high 
energy. We can conclude, there are very big dif- 
ferences in implementing the minijet components 
in the models. 

5. Properties of the models in the highest 
energy region 

In Fig.^T] the pseudorapidity distributions for 
charged hadrons according to DPMJET-II.3 are 
presented for energies between -^/s — 1 TeV and 
2000 TeV. The width of the distributions in- 
creases like the logarithm of the energy and also 
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Figure 7. Rapidity distribution of negatively 
charged hadrons in central S-S ans S-Ag colli- 
sions. The results of DPMJET-II.3 are compared 
with data from the NA-35 Collaboration fE3|. 



Figure 8. The jet transverse energy distri- 
bution is compared with data from the CDF- 
Collaboration Q|. The jets are found from the 
model events in the pseudorapidity region \r)\ < 
0.7 using a jet finding algorithm. 



the maximum of the curves rises like the loga- 
rithm of the energy. If we call the central re- 
gion around the two maxima the plateau, then we 
find the width of this plateau hardly to change 
with energy. Fig.[l^ presents the rise of the to- 
tal charged multiplicity with the cms energy ^fs 
according to DPMJET, QGSJET and SIBYLL. 
we find again, at low energies, where data are 
available, the models agree rather well. DPM- 
JET and SIBYLL agree in all the energy range 
shown. However, QGSJET above the energy of 
the TEVATRON extrapolates to higher energies 
in a completely different way. 

In Fig.|l3| we present for pp anf p-Air collisions 
the energy fractions K for B — B (baryon - an- 
tibaryon) and charged pion production. The cos- 
mic ray spectrum-weighted moments in p-A col- 
lisions are defined as moments of the F(xi a b) '■ 



Z i A — / {Xlab) 1 1 F[ A (xiab)dXlab 

Jo 



(3) 



mic ray energy spectrum and A represents both 
the target nucleus name and its mass number. 
In FigJlJ we present the spectrum weighted mo- 
ments for pion production in pp and p-Air colli- 
sions as function of the cms energy s/s per mi- 
cleon. We find all average values characterizing 
hadron production: the cross sections (Fig.|l|) , the 
average transverse momenta (FigjT^) the charged 



multiplicities (Fig 



1^), and the moments in Figs. 



Here —7 ~ -1.7 is the power of the integral cos- 



[14| and in Fig. [13 to change smoothly with en- 
ergy in most cases just like the logarithm of the 
energy. 

Comparison of the models after simulating 
the Cosmic Ray cascade 

First we present results of a comparison be- 
tween the cascade code HEMAS @ using DP- 
MJET as event generator and the cascade code 
CORSIKA II using VENUS as event generator 
Ga]. This comparison has been done for quan- 
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Figure 9. Comparison of transverse momentum 
cross sections according to DPMJET-II.3 at \fe 
= 0.63 TeV with collider data from the UA-1 
MIMI Collaboration |45| . The experimental date 
are represented by the parametrization given by 
the Experiment. 



Figure 10. Average transverse momenta of 
charged secondaries produced in pp and pp col- 
lisions calculated from DPMJET, QGSJET and 
SIBYLL (The latter two as given in Ref.g) as 
function of the center of mass energy y/s com- 
pared to date collected by the UA-1 Collabora- 
tion §§. 



tities of interest for the EAS-Top and MACRO 
experiments in the Gran Sasso Lab. The zenith 
angle is fixed at 31 degrees (MACRO/EAS-TOP 
coincidence direction). The e.m. shower size and 
muons above 1 TeV are sampled at 2000 meters 
a.s.l. (946 g/cm2 slant depth, 810 g/cm2 vert, 
depth). The calculations were done for primary 
protons, He nuclei and Fe nuclei with energies 
between 3 and 2000 TeV. Calculated are for each 
primary energy and particle (i) the e.m. shower 
profile, (ii)the Log(e.m. size) at EAS-TOP sam- 
pling depth (946 g/cm2), (iii)the distance muon- 
shower axis for E > 1 TeV muons, (iv) the muon 
decoherence for E > 1 TeV muons, (v) the number 
of muons per shower and (vi) the energy spec- 
trum of E > 1 TeV muons. In Figs. [l5| to |l8| 
we present two of these comparisons. A satis- 
factory agreement is found in these plots as well 
as in all other comparisons at different energies 
and with the other primary particles. Next we 



present two comparisons from the Karlsruhe code 
comparison j^] . The distributions choosen in this 
comparison are motivated by the interest of the 
KASKADEH experiment in Karlsruhe. In FigJlJ 
the Muon multiplicity distribution at ground level 
is calculated for primary protons of E = 10 15 eV. 
The calculation is done with the CORSIKA cas- 
cade code using 5 different event generators for 
the hadronic interactions. While again VENUS 
and DPMJET give distributions , which agree 
very well, it is found, that SIBYLL gives a very 
different distribution centered at smaller Muon 
number. 

In Figs. [2(] and ^l] ( The distributions were 
calculated using the CORSIKA shower code g 
with 5 different event generators for the hadronic 
interactions. ) Fe and p induced showers with 
energies of E = 10 14 and 10 15 eV are plotted 
in the log 10 -log 10 N e plane (Muon-number 
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Figure 11. The development of the pseudora- 
pidity distribution of charged hadrons produced 
in inelastic pp collisions in the in the center of 
mass energy range between y / s=l TeV and y/s = 
2000 TeV. 



-Electron-number plane). The distribution of 
events according to each of the 5 interaction mod- 
els for each energy and primary perticle is indi- 
cated by contours. Considering these plots calcu- 
lated with only one of the models, where Muon 
number is plotted over electron number, the im- 
pression is, that a simultaneous measurement of 
Muon-number and Electron number allows to de- 
termine the primary energy as well as the compo- 
sition of the primary component. In these plots 
we see, that for instance VENUS and DPMJET 
agree very well, but the contour according to 
SIBYLL for Fe projectiles of E — 10 15 eV over- 
laps the VENUS and DPMJET contours for p 
projectiles. From these differences between the 
models one can conclude, that at present the sys- 
tematic errors of the cascade calculations (and 
this are just the differences obtained using differ- 
ent models) prevent to identify safely the compo- 
sition of the primary component from such mea- 
surements. 



Figure 12. Rise of the charged multiplicity in in- 
elastic pp collisions according to DPMJET-II.3 in 
the center of mass energy range between y / s=0. 02 
TeV and <Js = 2000 TeV. At energies between 1 
and 100 TeV we plot also the average multiplici- 
ties according to SIBYLL and QGSJET as given 
in Rcf.g. 



6. Conclusions 

I would like to stress, more efforts are needed to 
extend the models used to simulate the hadronic 
interactions in the C.R. cascade up to the energies 
to be explored by the Auger Experiment. 

At least at collider energies, where data are 
available, these models should agree among them- 
selves and with the data. Disagreements to data 
like the ones seen in Fig.|| should be removed as 
soon as discovered. 

A much better understanding is needed how 
to calculate the minijet component. Certainly, 
the parton structure functions used for calculat- 
ing the minijet cross sections should correspond 
to the HERA measurements at small x. But this 
is certainly not the only problem. The differences 
in the extrapolation to higher energies of quanti- 
ties like average transverse momenta and charged 
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Figure 13. Laboratory energy fractions for B—B 
and pion production in pp and p-Air collisions 
according to DPMJET-II.3 as function of the 
(nucleon-nucleon) cms energy y^i. 



Figure 14. Spectrum weighted moments for pion 
production in pp and p-Air collisions as function 
of the (nucleon-nucleon) cms energy ^/s. 



multiplicities (see Figs. |l0| and ) in the three 
models implementing minijets are huge. These 
differences indicate, that much effort is needed to 
get a better understanding of the minijet compo- 
nent. 

Another question, where models disagree is the 
presence at high energy of an important soft com- 
ponent of hadron production like in the models 
with a supercritical pomeron. In minijet models 
all rise of the cross sections and of particle pro- 
duction at high energy is only due to the minijets. 

There are (even at energies, where collider data 
are available, see Fig. |l9|) large differences be- 
tween the models after simulating the C.R. cas- 
cade. We have to interpret these differences as 
the systematic errors of the cascade simulation. 
Such large differences could well prevent the in- 
terpretation of otherwise very interesting Cosmic 
Ray data. In future, C.R. results should always 
be interpreted using simulations with some differ- 
ent models. 

It might be dangerous, that at present many of 



the popular models are based on the same the- 
oretical foundations (and yet might differ very 
much in their results). To be on the safe side, it 
would be usefull to construct models based also 
on widely different theoretical concepts (for in- 
stance on the string fusion model |l9||). 

Finally, I would like to stress the need for new 
measurements of hadron production especially at 
the highest possible energies. In particular in the 
fragmentation region so important for the cosmic 
ray cascade, data (like Feynman-x distibutions) 
from the TEVATRON collider would be highly 
welcome. 
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Figure 19. Muon number distribution at ground 
level from proton induced showers with E = 10 15 
eV. The distributions were calculated using the 
CORSIKA shower code || with 5 different event 
generators for the hadronic interactions. 
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Projecting along the lines (m) and (E) one can 
estimate the energy and mass of the primary. 
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Figure 20. Contours in the log 10 -log 10 N e 
plane for p and Fe induced showers of E = 10 14 
and 10 15 eV . 
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